library(MatchIt)
library(cem)

##################################### Figure 7 ############## 
############################################################
###### Model 7: Using Matching
##matching by joint democracy

pre_match <- data %>% 
              select(challenger, target, year, ddyad, initiation, tgt_landneighbors2, backdown_terr1, backdown_nonterr1, 
                    joint_demo, defensetgt, offensechal, land_neighbor, tpopchal,tpoptgt,
                    pot_milicapchal, pot_milicaptgt, intra_wartgt, inter_war_totaltgt, 
                    solschangetgtdum, solschangechaldum, vicratio, capratio, tmid, tmid_sq,
                    tmid_cub)
       
pre_match <- na.omit(pre_match)
##run matching using Coarsened Exact Matching method
m.out <- matchit(joint_demo ~ defensetgt + log(tpopchal + 1) + log(tpoptgt + 1) +  offensechal +
                        intra_wartgt + inter_war_totaltgt,
                 data = pre_match, method = "cem")

summary(m.out)
plot(m.out)
###get matched data for model 
matched_data <- match.data(m.out)
names(matched_data)

mb1 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 + 
                      joint_demo +  defensetgt + offensechal +  land_neighbor + 
                      vicratio +  tmid +  tmid_sq + tmid_cub, 
                data = matched_data, weights = matched_data$weights,
                family = binomial(link = "logit"),
               prior.scale=2.5, prior.df=1)

#### Figiure  7(a)
coef_mb1 <- coef_clusterplot(ModelResults = mb1,
                            varnames = c("tgt_landneighbors2","backdown_terr1",
                                         "backdown_nonterr1", "joint_demo",
                                         "defensetgt", "offensechal",
                                         "land_neighbor","vicratio"),
                            data = matched_data, 
                            clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor", "Offensive alliance (challenger)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 

ggsave("./figures/coef_mb1.pdf", plot = coef_mb1,  width = 10, height = 8)

#### Model 8
## use negative pecae
rb2 <- bayesglm(initiation ~ negative_peacetgt +  backdown_terr1 +  backdown_nonterr1 + 
                       joint_demo +  defensetgt + offensechal +  land_neighbor + 
                       vicratio +  tmid +  tmid_sq + tmid_cub, 
                data = data, family = binomial(link = "logit"),
                prior.scale=2.5, prior.df=1)
## Figiure  7(b)
coef_rb2 <- coef_clusterplot(ModelResults = rb2,
                            varnames = c("negative_peacetgt","backdown_terr1",
                                         "backdown_nonterr1", "joint_demo",
                                         "defensetgt", "offensechal",
                                         "land_neighbor","vicratio"),
                            data = data, 
                            clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor", "Num. of states in negative peace",
                                 "Offensive alliance (challenger)",
                                 "Exp. victory prob. ratio")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 
ggsave("./figures/coef_rb2.pdf", plot = coef_rb2,  width = 10, height = 8)

## Medel 9 
##general backdown
rb3 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown1 + 
                       joint_demo +  defensetgt + offensechal +  land_neighbor + 
                       vicratio +  tmid +  tmid_sq + tmid_cub, 
                data = data, family = binomial(link = "logit"),
                prior.scale=2.5, prior.df=1)
## Figure 9(c)
coef_rb3 <- coef_clusterplot(ModelResults = rb3,
                             varnames = c("tgt_landneighbors2","backdown1",
                                          "joint_demo",
                                          "defensetgt", "offensechal",
                                          "land_neighbor","vicratio"),
                             data = data, 
                             clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (general)",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor",
                                 "Offensive alliance (challenger)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 

ggsave("./figures/coef_rb3.pdf", plot = coef_rb3,  width = 10, height = 8)


## Model 10
## both
rb4 <- bayesglm(initiation ~ negative_peacetgt +  backdown1 + 
                       joint_demo +  defensetgt + offensechal +  land_neighbor + 
                       vicratio +  tmid +  tmid_sq + tmid_cub, 
                data = data, family = binomial(link = "logit"),
                prior.scale=2.5, prior.df=1)
## Figure 9(d)
coef_rb4 <- coef_clusterplot(ModelResults = rb4,
                             varnames = c("negative_peacetgt","backdown1",
                                           "joint_demo",
                                          "defensetgt", "offensechal",
                                          "land_neighbor","vicratio"),
                             data = data, 
                             clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (general)",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor", "Num. of states in negative peace",
                                 "Offensive alliance (challenger)",
                                 "Exp. victory prob. ratio")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 
ggsave("./figures/coef_rb4.pdf", plot = coef_rb4,  width = 10, height = 8)

### Model 11: fixed effect logit model 
############## Figure 7(e) in the stata
# Model 11 was run using Stata. The stata code are provided below.

#      xtset ddyad year
#      *fixed effect
#      xtlogit initiation tgt_landneighbors2  backdown_terr1 backdown_nonterr1 ///
#       joint_demo defensetgt offensechal land_neighbor vicratio tmid tmid_sq tmid_cub if weaker_chal ==1, fe 	

#Save the stata results as dta file and read into R to make Figure 7(e)
library(readstata13)
fixlogit <- read.dta13("Model11.dta")
#Figure 7(e)
fixlogit <- fixlogit %>%
       dplyr::mutate(color_sig = ifelse(ci_lower < 0 & ci_upper > 0, "#ef8a62", "#2c7bb6"))
fixlogit <- fixlogit[1:8,]
fixedlogitstata <- ggplot(fixlogit, aes(x=var, y=coef, ymin=ci_lower, ymax=ci_upper)) +
       geom_pointrange(size=1, shape=16, colour=fixlogit$color_sig) +
       coord_flip() +
       geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
       xlab('') + ylab('') +  
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) +
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor", "Offensive alliance (challenger)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio")) 
ggsave("./figures/fixedlogitstata.pdf", plot =fixedlogitstata,  width = 10, height = 8)



### You can also use R for Model 11; The results are identical with the stata.
library(glmmML)
## fixed-effects logit 
logitfemod1 <- glmmboot(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 + 
                               joint_demo +  defensetgt + offensechal +  land_neighbor + 
                               vicratio +  tmid +  tmid_sq + tmid_cub,
                        family=binomial(link="logit"), data=data, cluster=ddyad)

## Figure 7(e): using R results
model = logitfemod1
modelcoef = model$coefficients
modelse = model$sd
ylo = modelcoef - 1.96* modelse
yhi = modelcoef + 1.96*modelse
names = names(model$coefficients)
dfplot = data.frame(names, modelcoef, modelse, ylo, yhi)
dfplot$names <- as.character(dfplot$names)

dfplot <- dfplot %>%
       dplyr::mutate(color_sig = ifelse(ylo < 0 & yhi > 0, "#ef8a62", "#2c7bb6"))
dfplot <- dfplot[1:8,]
fixedlogit <- ggplot(dfplot, aes(x=names, y=modelcoef, ymin=ylo, ymax=yhi)) +
       geom_pointrange(size=1, shape=16, colour=dfplot$color_sig) +
       coord_flip() +
       geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
       xlab('') + ylab('') +  
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) +
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor", "Offensive alliance (challenger)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio")) 
ggsave("./figures/fixedlogit.pdf", plot =fixedlogit,  width = 10, height = 8)


################################# Table 1 #################################
###########################################################################
#define function to get N
count <- function(x){
       length(which(!is.na(x)))
}
# all
summ <- function(x){ c(mean(x,na.rm = T),
                       sd(x, na.rm = T),
                       range(x,na.rm = T),
                       count(x)) }

##first at weaker dyad 
sumdf_weak <- data %>% 
       select(initiation, terrmidl, 
              tgt_landneighbors2, backdown_terr1, backdown_nonterr1, 
              joint_demo, defensetgt, offensechal, land_neighbor, 
              pot_milicapchal, pot_milicaptgt, intra_wartgt, inter_war_totaltgt, 
              solschangetgtdum, solschangechaldum, vicratio, capratio,
              negative_peacetgt, backdown1,tmid, tmid_sq, tmid_cub)

sumdf_weak <- sapply(sumdf_weak, function(x) summ(x))
sumdf_weak <- t(sumdf_weak)
sumdf_weak <- as.data.frame(sumdf_weak)
sumdf_weak$var <- rownames(sumdf_weak)
colnames(sumdf_weak)[1:5] <- c("mean", "sd", "min", "max", "N")
library(stargazer)


## Table 1
sumdf_weak <- as.data.frame(sumdf_weak)
stargazer(sumdf_weak, title = "Descriptive statistics", label = "des1",
          covariate.labels = c("Territorial claim initiation", "Territorial MID initiation",
                               "Num. of land neighbor (target)",
                               "Past yielding", "Past yielding (non-territory)",
                               "Joint democracy",
                               "Defensive alliance (target)",
                               "Offensive alliance (challenger)",
                               "Land neighbor",
                               "Alliance strength (challenger)",
                               "Alliance strength (target)",
                               "Civil war (target)",
                               "Interstate war (target)",
                               "SOLS change (target)",
                               "SOLS change (challenger)",
                               "Exp. victory prob. ratio",
                               "Capability ratio",
                               "Num. of states in negative peace",
                               "Past yielding (general)","Year since last initiation", 
                               "Year since last initiation(squared)",
                               "Year since last initiation(cub)"),
          type = "latex", digits = 1, out = "./figures/descrip1.tex")

## Figure A1: appendix (making a correlation plot)
cordf <- data %>% 
       select(initiation, tgt_landneighbors2, backdown_terr1, backdown_nonterr1, 
                         joint_demo, defensetgt, offensechal, land_neighbor, 
                         pot_milicapchal, pot_milicaptgt, intra_wartgt, inter_war_totaltgt, 
                         solschangetgtdum, solschangechaldum, vicratio, capratio)
colnames(cordf) <- seq(0, 15)

colnames(cordf)[1] <- "Y"
cor1 <- ggcorr(cordf,geom = "circle", nbreaks = 4, legend.position = "bottom", label = F,
               min_size = 0, max_size = 6) 
ggsave("./figures/cor1.pdf", plot = cor1,  width = 7, height = 7.2)



